%% Nyquist Frequency and Negative Frequencies


SampleRate=100; %Hz
Tvec=(0:1/SampleRate:1) %s

%We can play around with the frequency relative to the sampling rate

%RelFreq = 1/4
%RelFreq = 2/4
RelFreq = 3/4

Freq= RelFreq * SampleRate; %Hz

N=length(Tvec);
%This is our wave that will muliply the signal
Signal = exp(-i*2*pi*Freq*Tvec);

%This is how the above wave would look at our sampling frequency.  Note
%that RelFreq above 1/2 look like waves with relative frequency less than
%1/2. They can be seen to have frequencies of -(Freq-SampleRate)
figure
plot(Tvec,Signal,'o-')
xlabel('Time (s)')
ylabel('Signal')
hold on

% This is the actual wave at much higher sampleing rate.  Notice the where
%the red crosses the circles, which are the sampling points 
SampleRate=10000; %Hz
Tvec=(0:1/SampleRate:1) %s
Signal = exp(-i*2*pi*Freq*Tvec);
plot(Tvec,Signal)
legend('Sampled Wave','Actual Wave')


%% Resampling

%We can resample or recontruct the wave by padding the FFT since the sampling
% interval in time after inverse transform is 1/MaxFreq.  Any wave with frequency
%below the Nyquist frequency can be recontructed 

RelFreq = 1/4
SampleRate=100; %Hz
Freq= RelFreq * SampleRate %Hz

Tvec=(0:1/SampleRate:.5) %s

N=length(Tvec);
Signal = exp(2*pi*i*Freq*Tvec);

figure
plot(Tvec,Signal,'o-')

FS=fft(Signal);
%pad 
Npad=10;
FS=cat(2,FS,repmat(FS,[1 Npad])*0);
Tsub=1/SampleRate/(Npad+1)*(0:length(FS)-1);
Ssub=ifft(FS);
%The size and therfore normailzation is different on forward and inverse
%transforms, so we need to correct for that:
Ssub=Ssub*length(Ssub)/length(Signal);

hold on
plot(Tsub,real(Ssub),'-')
axis([.2 .3 -1 1])
xlabel('Time (s)')
ylabel('Signal')
legend('Sampled','Reconstructed')

%% Shifting FFT so Zero is at Center

close all

%Make a square pulse signal
N=100;
T=(0:N-1);
S=zeros(N,1);
S(T>30&T<50)=1;

figure
plot(T,S)
xlabel('T')
ylabel('Signal')
title('Signal')


%The Fourier Transoform
F=1/N*(0:N-1);
FTS=fft(S);
figure; hold on
plot(F,real(FTS))
plot(F,imag(FTS))
xlabel('F')
ylabel('FFT')
legend('real','imag')
title('Un-Shifted Transform')

% Shift zero to center
FTS_Shift=fftshift(FTS);

F_Shift=1/N*(-N/2:N/2-1)% For even number of points
%F_Shift=1/N*(-floor(N/2):floor(N/2))% For odd number of points

%The shifted Transform
figure
plot(F_Shift,real(FTS_Shift))
hold on
plot(F_Shift,imag(FTS_Shift))
xlabel('F')
ylabel('FFT')
legend('real','imag')
title('Shifted Transform')

%We could do some Fourier-domain filtering here.  Lets remove high frequencies

Ncut=20;
FTS_Shift(1:Ncut)=0;
FTS_Shift(end-Ncut+1:end)=0;

%The Filtered, shifted Transform
figure
plot(F_Shift,real(FTS_Shift))
hold on
plot(F_Shift,imag(FTS_Shift))
xlabel('F')
ylabel('FFT')
legend('real','imag')
title('Filtered, Shifted Transform')

%Inverse transform
FTS_Filter=ifftshift(FTS_Shift);
S_Filter=ifft(FTS_Filter);

%Plot our filtered signal
figure
plot(T,S_Filter)
xlabel('T')
ylabel('Signal')
title('Filtered Signal')

%notice Gibbs phenomenon
web('http://en.wikipedia.org/wiki/Gibbs_phenomenon')








